#!/bin/bash
#Quit on errors
set -e

# creat folder
echo "Create Output directory : out/"
mkdir -p out/

# Run this script by typing "q2_ad_script.sh | tee q2.lg"

# assembly
echo "1 assembly"
usearch11.0.667_i86linux64 -fastq_mergepairs *_1_1.fastq -reverse *_1_2.fastq -fastq_maxdiffs 10 \
          -fastq_pctid 20 -fastq_minmergelen 400 -fastq_maxmergelen 520 \
          -fastqout out/merged.fq -relabel @ -log out/merge.log

# Filter sequences
echo "2 Filter sequences"
usearch11.0.667_i86linux64 -fastq_filter out/merged.fq -fastq_maxee 1.0 -fastaout out/filtered.fa

# Find uniques
echo "3 Find uniques"
usearch10.0.240_i86linux32 -fastx_uniques out/filtered.fa -fastaout out/uniques.fa \
          -relabel Uniq -sizeout

# Cluster OTUs
echo "4 Cluster OTUs"
usearch11.0.667_i86linux64 -cluster_otus out/uniques.fa -minsize 2 -otus out/otus.fa \
          -relabel Otu

# OTU table
echo "5 OTU table"
usearch10.0.240_i86linux32 -otutab out/merged.fq -otus out/otus.fa \
          -notmatched out/notmatched.fa -otutabout out/otutab.txt

# Normalize OTU table
echo "6 Normalize OTU table"
usearch10.0.240_i86linux32 -otutab_norm out/otutab.txt -sample_size 50000 \
          -output out/otutab_norm.txt

#Uncross
echo "7 Uncross"
usearch10.0.240_i86linux32 -uncross out/otutab_norm.txt -tabbedout out/out.txt \
          -report out/rep.txt -otutabout out/result.txt

# Trim OTU table
echo "8 Trim OTU table"
usearch10.0.240_i86linux32 -otutab_trim out/result.txt -min_otu_freq 0.001 \
          -output out/trimmed.txt

col_otu=$(awk '{print NF; exit}' out/otutab_norm.txt)
col_otu_trim=$(awk '{print NF; exit}' out/trimmed.txt)
if [ "$col_otu" != "$col_otu_trim" ]; then
	echo "$(tput bold)$(tput setaf 1)#####################################################################################"
	echo "                         WARNING: Different number of samples !!!!!                          "
	echo "#####################################################################################$(tput sgr0)"
	exit 0
fi
 
echo $'\n'
echo "$(tput setaf 3)#######################################################################"
echo '08: Grep'
echo "#######################################################################$(tput sgr0)"
awk '{print $1 }' out/trimmed.txt > out/otus2trim.txt
usearch10.0.240_i86linux32 -fastx_getseqs out/otus.fa -labels out/otus2trim.txt \
          -fastaout out/otus_trimmed.fa

echo $'\n'
echo "$(tput setaf 3)#######################################################################"
echo '09: Import OTUs (fna) to Qiime2'
echo "#######################################################################$(tput sgr0)"
source /Users/local/bin/bioconda/activate qiime2-2020.2
qiime tools import --input-path out/otus_trimmed2.fa \
                   --output-path out/otus_trimmed2.qza \
                   --type 'FeatureData[Sequence]'

echo $'\n'
echo "$(tput setaf 3)#######################################################################"
echo '10: Assign taxa'
echo "#######################################################################$(tput sgr0)"
time qiime feature-classifier classify-consensus-blast \
           --i-query out/otus_trimmed2.qza \
           --i-reference-reads /home/adly/qiime/qiime2_2021/99_otus_16S.qza \
           --i-reference-taxonomy /home/adly/qiime/qiime2_2021/majority_taxonomy_all_levels.qza \
           --o-classification out/taxonomy_blast_silva2.qza \
           --p-perc-identity 0.97 \
           --p-min-consensus 0.8 \
           --p-maxaccepts 1

echo $'\n'
echo "$(tput setaf 3)######################################################################"
echo '13: Align Sequences'
echo "######################################################################$(tput sgr0)"
qiime alignment mafft --i-sequences out/otus_trimmed2.qza \
                --p-n-threads 64 \
                --o-alignment out/otus_trimmed_aligned2

echo $'\n'
echo "$(tput setaf 3)######################################################################"
echo '14: Phylogeny tree'
echo "######################################################################$(tput sgr0)"
qiime phylogeny fasttree --i-alignment out/otus_trimmed_aligned2.qza \
                --p-n-threads 64 \
                --o-tree out/Tree2

echo $'\n'
echo "$(tput setaf 3)######################################################################"
echo '15: midpoint tree'
echo "######################################################################$(tput sgr0)"
qiime phylogeny midpoint-root --i-tree out/Tree2.qza \
      --o-rooted-tree out/Tree_midpoint2

      
qiime tools export --input-path out/Tree_midpoint2.qza --output-path out/Tree_midpoint2
qiime tools export --input-path out/taxonomy_blast_silva2.qza --output-path out/taxonomy2

#calculate number of reads per samples
usearch11.0.667_i86linux64 -alpha_div trimmed.txt -output alpha.txt -metrics reads

#No of OTUs persample and richness
usearch11.0.667_i86linux64 -alpha_div trimmed.txt -output alpha_richness.txt -metrics richness

#convert the trimmed file to biom
biom convert -i out/trimmed.txt -o out/trimmed_json.biom --table-type="OTU table" --to-json

# Import biom file to qiime2
qiime tools import   --input-path out/trimmed_json.biom   --type 'FeatureTable[Frequency]'   --input-format BIOMV100Format   --output-path out/feature-table-2.qza

# calculate pielou_e_evances
qiime diversity alpha   --i-table feature-table-2.qza   --p-metric pielou_e  --o-alpha-diversity R1_R2-pielou_e-observed_otus_vector.qza

qiime tools export --input-path R1_R2-pielou_e-observed_otus_vector.qza --output-path R1_R2-pielou_e-observed_otus_vector2.qzaR1_R2-pielou_e-observed_otus_vector.tsv


# calculate Shannon
qiime diversity alpha   --i-table feature-table-2.qza   --p-metric shannon   --o-alpha-diversity R1_R2-shannon-observed_otus_vector.qza

qiime tools export --input-path R1_R2-shannon-observed_otus_vector.qza --output-path R1_R2-shannon-observed_otus_vector.tsv

# Calculate Simpson
qiime diversity alpha   --i-table feature-table-2.qza   --p-metric simpson   --o-alpha-diversity R1_R2-simpson-observed_otus_vector.qza

qiime tools export --input-path R1_R2-simpson-observed_otus_vector.qza --output-path R1_R2-simpson-observed_otus_vector.tsv


echo "###################################################################################"
echo "###################################################################################"
echo "                                      END SCRIPT                                   "
echo "###################################################################################"
echo "###################################################################################"

